-
Notifications
You must be signed in to change notification settings - Fork 2
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Use time-of-flight workflow in Dream data reduction #125
Conversation
"metadata": {}, | ||
"outputs": [], | ||
"source": [ | ||
"grouped_dspacing = workflow.compute(IofDspacingTwoTheta)\n", | ||
"grouped_dspacing = workflow.compute(IofDspacingTwoTheta).bins.concat(\"event_time_zero\")\n", |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Does this mean that the result has a new dimension?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, we expect NeXus data to have the event_time_zero
dimension, because it records events from more than one pulse.
The data we have in the geant4/mcstas files isn't very realistic, it has values for tof that exceed 71ms.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
But shouldn't the time have been removed at some point during conversion to d-spacing?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is gone after the latest update. We no longer bin the data in event_time_zero
, we keep it as an event coord.
# Add a event_time_zero coord for each bin, but not as bin edges, | ||
# as all events in the same pulse have the same event_time_zero, hence the `[:2]` | ||
da.coords["event_time_zero"] = ( | ||
sc.scalar(1730450434078980000, unit="ns").to(unit="us") + da.coords["tof"] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What is this number?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's a random date in the past. The actual value does not matter.
I added a comment.
src/ess/dream/io/geant4.py
Outdated
|
||
period = (1.0 / sc.scalar(14.0, unit="Hz")).to(unit="us") | ||
# Bin the data into bins with a 71ms period | ||
da = da.bin(tof=sc.arange("tof", 3) * period) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why 2 bins? Does it just so happen that the simulation contains 2 pulses?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It is 2 bins because the tofs spill over into the next pulse, but not in the one after that.
But in principle, you could have that in another file.
I changed to compute a npulses
based on the max tof and the period.
src/ess/powder/conversion.py
Outdated
def compute_detector_time_of_flight( | ||
detector_data: DetectorData[RunType], tof_workflow: TofWorkflow | ||
) -> TofData[RunType]: | ||
wf = tof_workflow.pipeline.copy() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
How expensive is this? Does it copy the cached intermediate results?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I always assumed it made a cheap shallow copy?
I think it's using https://networkx.org/documentation/stable/reference/classes/generated/networkx.Graph.copy.html?
src/ess/powder/conversion.py
Outdated
def set_monitor_zeros_to_nan( | ||
monitor: TofMonitorData[RunType, MonitorType], | ||
) -> TofMonitorDataZerosToNan[RunType, MonitorType]: | ||
inds = monitor.values == 0.0 | ||
monitor.values[inds] = np.nan | ||
return TofMonitorDataZerosToNan[RunType, MonitorType](monitor) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looking back at the discussion about granularity, can this be merged with the above provider? Do we need to expose the data without NaNs here? Or should it be part of convert_monitor_to_wavelength
?
Also, this provider modifies its input which can lead to surprises down the line. Avoiding this would be easiest by merging it into compute_monitor_time_of_flight
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah, I changed my mind at least twice on this one.
In the end, I made it into separate steps because I thought maybe not all workflows or instruments would want to NaN the zero counts, and it would make it easier to insert a different provider.
But I agree about making the graph more messy and modifying the data in place.
src/ess/snspowder/powgen/workflow.py
Outdated
detector_data: | ||
Data from the detector. | ||
""" | ||
return TofData[RunType](detector_data) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we avoid this dummy provider by changing extract_raw_data
to return TofData
instead of DetectorData
?
src/ess/dream/io/geant4.py
Outdated
npulses = int((da.bins.coords["tof"].max() / period).value) | ||
da = da.bin(tof=sc.arange("tof", npulses + 1) * period) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
npulses = int((da.bins.coords["tof"].max() / period).value) | |
da = da.bin(tof=sc.arange("tof", npulses + 1) * period) | |
npulses = int((da.bins.coords["tof"].max() / period).ceil().value) | |
da = da.bin(tof=sc.arange("tof", npulses) * period) |
So that npulses
actually is the number of pulses. Without this, I think [:npulses]
below is incorrect.
@jl-wynen see update. It's now cleaner after the latest changes in The CI will fail until we release |
I also updated the use of metadata utilities after they were updated in scippneutron. EDIT: Hmm it seems I should have waited for #124 ... |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this notebook supposed to be shown in the docs? If so it needs to be added to the ToC and needs some descriptions. If not, please move it somewhere else, e.g., /tools
because it gets executed when building the docs.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good point, I think I will actually move this to the essreduce
docs, that shows in a generic way how to create a lookup table. I think moving it to a tools
folder is good for now.
"wf[time_of_flight.LookupTableRelativeErrorThreshold] = 0.02\n", | ||
"# wf[time_of_flight.PulsePeriod] = 1.0 / sc.scalar(14.0, unit=\"Hz\")\n", | ||
"# wf[time_of_flight.PulseStride] = 1\n", | ||
"# wf[time_of_flight.PulseStrideOffset] = None" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"wf[time_of_flight.LookupTableRelativeErrorThreshold] = 0.02\n", | |
"# wf[time_of_flight.PulsePeriod] = 1.0 / sc.scalar(14.0, unit=\"Hz\")\n", | |
"# wf[time_of_flight.PulseStride] = 1\n", | |
"# wf[time_of_flight.PulseStrideOffset] = None" | |
"wf[time_of_flight.LookupTableRelativeErrorThreshold] = 0.02" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I uncommented them instead.
src/ess/dream/data.py
Outdated
|
||
|
||
def tof_lookup_table_high_flux() -> str: | ||
"""Path to a HDF5 file containing a lookup table for high-flux ToF.""" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
How was this created? using McStas or ToF? (Please explain in the docstring)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done.
Co-authored-by: Jan-Lukas Wynen <[email protected]>
In this PR, we incorporate the time-of-flight computation workflow into the DREAM data reduction workflow.
This results in much improved results where peaks in d-spacing are narrow and symmetrical, and lines are vertical along the 2theta dimension.
Before:


After:

